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Abstract. In this article, we motivate the detailed comparison of the physical prop- 
erties of individual configurations of a ferroelectric solid solution as a means toward 
developing first principles models for these systems. We compare energies, dielectric 
constants eoo, mode effective charges of local polar distortions, and the zero tempera- 
ture piezoelectric behavior of several ordered PbsGeTe4 supercells. Cluster expansions 
of these properties show the importance of second-neighbor effects, which can be re- 
lated to symmetry-preserving relaxation and its effect on the symmetry breaking polar 
instabilities. 



INTRODUCTION 

Ferroelectric solid solutions are of great technological importance. For exam- 
ple, the largest piezoelectric response are found in mixed ferroelectrics, such as 
Pb(Zri_ a; Ti a ;)03. Recently, giant piezoelectricity was discovered in the relaxor fer- 
roelectric systems Pb(A 1 / 3 Nb2/3)0 3 -PbTi03, (A = Zn, Mg). [1] Understanding the 
physics of ferroelectric solid solutions on the microscopic level would be of great 
theoretical interest and could also point to new ways to tune their piezoelectric 
response and other properties. 

Ab initio calculations have proved successful in relating properties of stoichio- 
metric ferroelectrics to phenomena on the atomic level [2]. Structural parameters 
[2,3], dielectric constants, effective charges, phonon dispersion relations [4,5], and 
polarizations [6] have been calculated from first principles. To predict the behavior 
at finite temperature, first-principles models have been constructed. [7—11]. Gener- 
ally, these models are based on a vector representation of the local polar distortions 
responsible for the ferroelectric phase transition. An electric dipole moment is as- 
sociated with each local distortion. At long distance, the interaction between local 
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distortions is dipole-dipole; significant corrections appear at short range. Anhar- 
monic terms, elastic constants and strain coupling to local distortion also appear, 
determining the ground state and the nature of the phase transitions in the mod- 
els. [12] The models obtained allow for simulation of ferroelectric phase transitions, 
where the Curie point is generally in good agreement with experiment. The models 
also allow piezoelectricity and dielectric functions to be computed [13,14]. In the 
following, we discuss the modeling of ferroelectric and piezoelectric behavior in solid 
solutions, and lay part of the groundwork for a Pb^^Ge^Te model by exploring the 
configuration dependence of certain quantities essential for constructing the model. 

WHY COMPARE CONFIGURATIONS? 

Consider a ferroelectric solid solution. The physical properties of this solid solu- 
tion are the ensemble averages of the properties of the individual ensemble mem- 
bers. Thus one needs to be able to compute these physical properties for individual 
ensemble members. A given ensemble member (configuration) can be treated as 
if it were an actual crystal structure and modeled in the same manner as for a 
stoichiometric ferroelectric. 

While models for a few chosen small unit cell configurations are obtained using 
this approach, some model for the general configuration is necessary in order to 
properly take ensemble averages. Given the nature of the model described in the 
Introduction, we expect that in a model for the general configuration, there will 
again be polar instabilities, a local basis for these distortions, and interactions 
between the local distortions. However, now the local polar distortions themselves 
and their interactions will be site-dependent. In principle, the site-dependence of 
any quantity can be described by a cluster expansion. The models obtained for the 
series of individual small unit cell configurations can be used to obtain the cluster 
expansions for a model that is valid for all configurations. 

In modeling stoichiometric ferroelectrics, two kinds of truncation in the mod- 
els are necessary to prevent an explosion of terms: the order of the expansion of 
the energy in powers of the local distortion variable and the range of the local 
interactions. The cluster expansion approach necessitates a third kind of physi- 
cally motivated truncation: the range of the cluster expansion. Where possible, 
terms should be found whose configuration dependence is unimportant and then 
kept constant for all configurations. For those terms where a cluster expansion is 
necessary, it should be truncated at the shortest possible range. Toward this end, 
it is thus very important to establish and understand the configuration dependence 
of those quantities that determine the model parameters. 

Pbi-zGe^Te provides an excellent prototype system for investigating the form 
of first-principles models in solid solutions using the approach just outlined. For 
all compositions above a critical composition x w 0.005 [15], Pb 1 _ a; Ge ;E Te under- 
goes a transition from a cubic phase for T > T c (x) to a rhombohedral phase at 
T < T c (x) [16]. Because the endmembers of the solid solution series have the 



rocksalt structure with only 2 atoms per primitive cell, a large number of differ- 
ent mixed configurations can be investigated without the need for excessively large 
supercells. To simplify the problem further, we focus on configuration dependence 
at a fixed composition, Pbo.75Geo.25Te. This composition was chosen to be near 
the low Ge concentration regime of physical interest, while allowing for a variety 
of supercells with only 8 atoms. We present results on the five 8-atom supercells 
of highest symmetry, shown in Figure 1. In each configuration studied, all Ge ions 
are translationally equivalent and therefore have the same environment, greatly 
simplifying the analysis. 



FIRST PRINCIPLES METHODS 



In our study of the structural phase transition in the cP8 configuration [10], 
we give full details of the first principles calculations. Briefly, we performed all 
calculations using density functional theory within the local density approximation 
(LDA). Ab initio pseudopotentials were used for the ions and a plane wave basis set 
with cutoff energy 300 eV was used for the Kohn-Sham eigenfunctions. Total en- 
ergies were calculated via conjugate gradients optimization using the CASTEP 2.1 
package [17]. The CASTEP 2.1 package was modified [18] to do variational linear 
response [5] calculations of force constants, Born effective charges and dielectric 
constants. Brillouin zone averages for the cP8 configuration were performed using 
a 4 x 4 x 4 Monkhorst-Pack set. For the other configurations, the same k point 
grid was used, folded into the corresponding Brillouin zone. 




FIGURE 1. Pb3GeTe4 configurations investigated. Each has 8 atoms per primitive unit cell. 
(The Pearson symbols specify the number of atoms in the conventional unit cell.) 



SYMMETRY PRESERVING RELAXATION 



We found two distinct ways in which the total energies of the configurations 
that we studied could be lowered with respect to the energy of the structures 
with all atoms fixed on rocksalt positions: symmetry-preserving relaxation and 
symmetry-breaking polar instabilities. The existence of symmetry-preserving re- 
laxation follows from group theory. Except for the cP8 configuration, all of the 
configurations have one or more identity irreducible representations (irreps) among 
its normal modes. An identity irrep implies an energy term that is linear in the 
corresponding mode and thus a lowering of the energy by relaxation in the ionic 
displacement subspace spanned by these modes. This relaxation neither lowers the 
symmetry nor causes net polarization in the configurations studied. The other nor- 
mal modes in the systems have non-identity irreps; thus to lowest order, the total 
energy is quadratic in these modes. The existence of instabilities (normal modes 
whose harmonic coefficient is negative) with their concurrent symmetry breaking 
is a system-specific phenomenon, as is the polar nature of these instabilities. This 
section focuses on symmetry-preserving relaxation; the next section focuses on 
symmetry-breaking polar instabilities. 

Before investigating symmetry-preserving relaxation, we relaxed each configu- 
ration of Figure 1 with respect to strain, holding all ions fixed at ideal rocksalt 
positions. At this point, the total energies all differed by less than 10 meV/ (8 atom 
cell) and the (pseudocubic) lattice parameters all differed by less than 0.06%. Had 
all subsequent calculations been performed with respect to the same pseudocubic 
lattice parameters, e. g. the cP8 value of a = 6.275 A, the strain energies involved 
would been negligible. 

Each configuration was allowed to undergo symmetry preserving relaxation, with 
its strain held fixed. The relaxed atomic coordinates can be simply described (to 
within 0.03 A) by displacement of the Te ions in each Pb-Te-Ge chain segment 
0.12 A toward to the Ge ion, with all other ions held fixed. [19] Qualitatively, re- 
laxation in Pbi-^Ge^Te can be regarded as the result of the size mismatch between 
the smaller Ge and the larger Pb ions. The energy of the cP8 configuration, for 
which no symmetry-preserving relaxation occurs, was 171 meV/cell higher than for 
the relaxed oC16 configuration, in which every Te is displaced. 

The energies of the relaxed configurations were fit to a pair-only cluster expan- 
sion out to fourth cation-cation neighbors. In meV/cell, the total energy of the 
configurations tested is given by: 

E = E - 2.97?! + 2QAN 2 - 1.5 N 3 + 1.37V 4 , (1) 

where iVj is the average number of i'th neighbor Ge ions per Ge ion. The 
cation sublattice in Pbi_ x GeTe is fee; each cation has 12 first neighbors cations 
at (a/2, a/2, 0), etc., 6 second neighbors at (0,0, a), etc., 24 third neighbors at 
(a/2, a/2, a), etc., and 12 fourth neighbors at (a, a, 0), etc., where a is the conven- 
tional fee lattice parameter. The sensitivity of energy to N 2 is an order of magnitude 



larger than for the other terms, illustrating the importance of the above-mentioned 
relaxation. The Ge and the Pb in a linear Ge-Te-Pb segment where the Te under- 
goes significant relaxation are second neighbors; the more Ge-Te-Pb segments, the 
fewer Ge-Te-Ge segments, the lower N2, and the lower the relaxed energy. 

The expression 1 is a fit of 5 relaxed cell energies to five parameters. It is 
also important to get an estimate of the predictive power of a cluster expansion 
truncated at a smaller number of terms. We did a series of cluster expansions in 
terms of N\, N2 and N% to each subset of 4 configurations [20] in Figure 1, and used 
the resulting expression to calculate the energy of the configuration that was left 
out. The rms error was 16 meV/cell and the maximum error was 20 meV/cell. This 
compares with the root mean square variance of 57 meV/ cell for the relaxed energies 
of the set of configurations studied, showing modest but not excellent predictability 
for a cluster expansion out to third neighbor. In what follows, we present results 
based on three or fewer configurations and the cluster expansions are truncated 
where the number of parameters matches the number of unknowns. Calculations 
on additional configurations would be necessary to make any meaningful test of 
predictability. 

SYMMETRY BREAKING INSTABILITIES 

Next, we looked at energy lowering via symmetry-breaking distortions. For each 
configuration, we then calculated the normal modes at the zone center. We per- 
formed linear response force constant calculations at q = for each relaxed config- 
uration and then diagonalized the corresponding dynamical matrix to obtain the 
normal modes. For each configuration, there were symmetry-breaking polar insta- 
bilities, with a three- dimensional vector representation centered on the Ge ion. [21] 
Further q 7^ calculations on cP8 and tP8 showed that instabilities throughout the 
Brillouin zone could be well- represented by a local (lattice Wannier function [24]) 
basis with one vector on each Ge ion. This supports the Ge off-centering model for 
ferroelectricity in Pbi-^Ge^Te [22] and gives strong evidence that the correct form 
for a first-principles model for arbitrary configurations will have one vector per Ge 
ion. 

In order to appropriately model the long-range physics within our models, both 
the electronic dielectric tensor and the polarization associated with the local 
instabilities must be determined. These are linear response functions of the relaxed 
high-symmetry structures. We have calculated the dielectric tensors and the 
Born effective charges for the cP8, tP8, and tI16 configurations. The Born effective 
charge tensors Z* and the normal mode ionic displacement pattern Ujp associated 
with the zone-center instability which transforms like the vector component a were 
used to calculate the so-called mode effective charge tensor through 

Ki = E(3W«tf)7 ( 2 ) 

31 



The normal modes were all normalized such that J2j \ u ji3\ 2 was 1 ^- 2 - 

The results for the and ~Z* tensors are shown in Table 1. Fitting to a pairwise 
cluster expansion out to second neighbor, we obtain: 

( eoo ) Qa = 42.5 - 0.1757V!,, - 0.1N 1± + 1.7 N n + 0.2N 2± , (3) 

where Ni±_ is the mean number of i'th neighbor Ge ions per Ge ion in a direction 
perpendicular to a and JVjii is the mean number of i'th neighbor Ge ions per Ge 
ion in a direction having nonzero component along a. Likewise, the mode effective 
charge is given (in units of eA) by 

T aa = 17.31 - 0.147?!!! - 0.127V 1± - 2.17JV 2 || - 0.22iV 2± . (4) 

The relative configuration dependence of Z is larger than that of e^. The three 
configurations all have diagonal dielectric and mode effective charge tensors; it 
should be noted that off-diagonal terms will be nonzero in the general (asymmetric) 
configuration. For the mode effective charge tensors, the more Ge-Te-Pb or Pb- 
Te-Ge segments along z, the lower N 2 \\ and the higher Z zz . We have previously 
shown how this increase can be explained in terms of local bonding [23]. The 
differences in the Born effective charges of individual ions between configurations 
is relatively small; it is primarily the difference in the normal mode eigenvector 
between configurations that leads to the increase in Z . In particular, in Ge-Te-Pb 
chain segments, relaxation of the Te away from the Pb leads effectively to an off- 
centering instability of the Pb ion that does not occur when the Pb-Te distance is 
smaller. [23] The active participation of the Pb ion helps to increase Z . 

TABLE 1. Electronic dielectric constant and mode effective charges in three 
Pb 3 GcTc4 configurations. 



Component 


cP8 


tP8 


tI16 


(foo)xa; 


46.7 


45.6 


46.3 


( e oo)yy 


46.7 


45.6 


46.3 


(^oo)zz 


46.7 


42.9 


43.3 


XX 


12.10 


11.97 


12.54 




12.10 


11.97 


12.54 




12.10 


15.95 


16.44 



PIEZOELECTRICITY AT ZERO TEMPERATURE 

The above calculations were for relaxed high symmetry states. We now turn our 
attention to properties of the fully distorted LDA ground states. We have previ- 
ously reported the piezoelectric tensors at zero temperature for the cP8 and tP8 



configurations [23]. In this section, we give further results for the tI16 configura- 
tion, write a simple cluster expansion for the zero-temperature piezoelectricity in 
Pbo.75Geo.25Te, and use this expansion to estimate the piezoelectric tensor of the 
disordered system. 

TABLE 2. LDA ground state of Pb3GeTe4-tI16. Monoclinic, space group Cm, a 
= 8.984 A, b = 8.876 A, c = 7.738 A, (3 = 54.08°. 



Atom 


Wyckoff position 


X 


y 


z 


Pb 


2(a) 


0.5034 





0.0061 


Pb 


4(b) 


0.5020 


0.2489 


0.5000 


Ge 


2(a) 


0.0291 





0.0069 


Te 


2(a) 


0.2355 





0.5158 


Te 


2(a) 


0.7578 





0.4753 


Te 


4(b) 


0.2391 


0.2383 


0.9976 



The LDA ground state for the tI16 configuration is given in Table 2. In Table 3, 
we give the full piezoelectric tensors for the LDA ground states of the cP8, tP8 
and tI16 configurations. In each case, the orientation is with respect to the axes of 
Figure 1, and the symmetry has been broken such that the Ge displacement and 
the polarization lie in the +{111} quadrant. 

TABLE 3. Comparison of piezoelectric strain tensors for the ground states of 
three Pb3GeTe4 configurations (in C/m 2 ). Components in parentheses are equal 
to other components via symmetry. Components which do not appear in the table 
are related by symmetry to those that do; e.g. e 2 2 = en for each configuration 



Component 


cP8 


tP8 


tI16 


en 


1.8 


2.5 


2.4 


ei2 


-0.6 


-1.2 


-1.2 


ei3 


(-0.6) 


-0.9 


-0.6 


e u 


-0.2 


-0.6 


-0.5 


ei5 


1.0 


1.2 


0.2 


ei6 


(1.0) 


1.1 


1.1 


e 3 i 


(-0.6) 


-0.3 


-0.5 


e33 


(1.8) 


5.1 


6.5 


e 3 4 


(1.0) 


2.2 


6.4 


e36 


(-0.2) 


-0.5 


-0.7 



If only pair terms are included in the cluster expansion for piezoelectricity, then 
the piezoelectric tensor is given by 

e = e + ^iV d e d , (5) 

d 



where e is the piezoelectric tensor, e is a configuration-independent term, {d} is 
the set of cation-cation separations, N d is the average number of Ge neighbors per 
Ge ion at separation d, and is the correction to the piezoelectric tensor due to 
neighbors at separation d. We are now dealing with a property of Pbo.75Geo.25Te at 
zero temperature; the discussion of piezoelectric tensor in this section only applies 
to the subensemble of zero temperature structures that have been poled by a field 
in the +111 direction. The subensemble has rhombohedral symmetry; the average 
piezoelectric tensor has only four independent components: en, eu, and e^. 

Symmetry constrains the form of the individual ed. Consider for example the 
effect of adding one Ge-Ge pair separated by (a/2, a/2,0) to a perfectly rhombo- 
hedral configuration. The global symmetry is broken to monoclinic. By writing 
the appropriate symmetrized form of the piezoelectric tensor for a monoclinic sys- 
tem and subtracting that for a rhombohedral system, the correct symmetry for 
e (a/2,a/2,o) is obtained. It is straightforward to show that the tensors ed for two dif- 
ferent d are related if and only if the two d are equivalent by symmetry under the 
rhombohedral group. In the cation sublattice, the 12 first neighbors break under 
rhombohedral distortion into 2 groups of 6, while the 6 second neighbors are all 
equivalent. 

The cP8, tP8 and tI16 configurations are sufficient to separate the first neighbor 
and second neighbor contributions to the piezoelectric tensor, but insufficient for 
determining the independent contributions of the two kinds of first neighbor. In 
Table 4, we give e and the contributions e( / 2 , a /2,o) + e( a /2 -a/2,0) and e (0 , ,a)- Both 
the tensors e( /2, /2,o) + e (a/2 ,-0/2,0) an d e (o,o,a) have the same symmetry as that for 
the tP8 and tI16 ground state piezoelectric tensors. 

The dominant terms reflect the piezoelectric components that change the most 
from configuration to configuration. The lack of second neighbor Ge pairs along 
(0,0, a) tends to increase the value of e 33 , while the lack of Ge pairs separated 
by (±a/2, ±a/2, 0) and (0,0, a) both increase the value of e 34 . The physics of the 
(0, 0, a) pairs is clear: the presence of a Ge-Te-Pb chain segment means that there is 
relaxation of the Te atom joining them and thus a weakening of the instability that 
transforms like z. A weak instability implies large response [23,14]. The e 33 and 
e 34 components are the ones enhanced because the instability along z is effectively 
coupled most strongly to these components. [10]. The exact source of the large 
influence of first neighbor pairs on e 34 has yet to be determined. 

Given the data in Table 4, it is possible to estimate the piezoelectric tensor of 
disordered Pbo.75Geo.25Te. We assume for present purposes that all members of 
the poled subensemble are equiprobable. Then there are no spatial correlations 
and iVd = 0.25 for all d. Using the values in Table 4 and Eq. 5, we estimate for 
disordered Pbo.75Geo.25Te that en = 5.9, ei 2 = —0.9, eu = —1.0, and e±5 = 3.7, 
reproducing the tensor form expected for rhombohedral symmetry. The estimated 
values of the components suggest that the piezoelectric response of the cP8 con- 
figuration is particularly unrepresentative for the disordered system. Because of 
the strong dependence of piezoelectric response on the magnitude of instability 
in a system and the configuration dependence of the magnitude of local instabil- 



TABLE 4. Terms in cluster expansion of Pbo.75Geo.25Te piezoelectric tensor (in 
Cjm 2 ) 



Component 


e 


e (a/2,a/2,0) + e (a/2,-a/2,0) 


e (0,0,a) 


en 


7.7 


0.1 


-0.3 


ei2 


-1.0 


0.2 


0.3 


ei3 


(-1.0) 


-0.1 


0.0 


ei4 


-1.2 


0.0 


0.1 


eis 


5.6 


0.5 


0.4 


ei6 


(5.6) 


0.0 


0.0 


e3i 


(-1.0) 


0.1 


-0.1 


e33 


(7.7) 


-0.7 


-2.4 


e34 


(5.6) 


-2.2 


-2.7 


e36 


(-1.2) 


0.1 


0.2 



ity, piezoelectric results on single high symmetry supercells of mixed ferroelectrics 
should be treated with caution. 

CONCLUSIONS 

This work describes the necessity of cluster expansions in first principles mod- 
els for piezoelectricity and ferroelectricity in solid solutions and the importance 
of comparing the properties of different configurations. In the specific example of 
Pbo.75Geo.25Te, strong configuration dependence is found for relaxed energies, the 
nature of the local polar instabilities, and the polarization associated with these 
instabilities. A good model for ferroelectricity in a mixed system must be able to 
account for this variability. We show how results of the zero temperature piezo- 
electricity on several ordered supercells can be used to obtain an estimate of the 
zero temperature piezoelectricity of a completely disordered cell. 

ACKNOWLEDGEMENTS 

This work was supported by ONR N00014-97-J-0047. 

REFERENCES 

1. S.-E. Park and T. R. Strout, in 1996 IEEE Ultrasonics Symposium Proceedings 
(New York: IEEE, 1996), v.2, p. 935; J. Appl. Phys. 82, 1804 (1997), and references 
therein. 

2. R. E. Cohen, Nature 358, 136 (1992). 

3. D. J. Singh, Phys. Rev. B 52, 12559 (1995). 



4. S. Baroni, P. Giannozzi and A. Testa, Phys. Rev. Lett. 58, 1861 (1987). 

5. X. Gonze, D. C. Allan and M. P. Teter, Phys. Rev. Lett. 68, 3603 (1992). 

6. R. D. King-Smith and D. Vanderbilt, Phys. Rev. B 47, 1651 (1993). 

7. K. M. Rabe and J. D. Joannopoulos, Phys. Rev. Lett. 59, 570 (1987). 

8. W. Zhong, D. Vanderbilt, and K. M. Rabe, Phys. Rev. Lett. 73, 1861 (1994); Phys. 
Rev. B 52, 6301 (1995). 

9. K. M. Rabe and U. V. Waghmare, Ferroelectrics 164, 15 (1995); U. V. Waghmare 
and K. M. Rabe, Phys. Rev. B 55, 6161 (1997). 

10. E. Cockayne and K. M. Rabe, Phys. Rev. B 56, 7947 (1997). 

11. H. Krakauer, R. Yu, C.-Z. Wang, K. M. Rabe and U. V. Waghmare, unpublished 
(cond-mat/9710088); H, Krakauer, R. Yu, C. Z. Wang, and C. Lasota, Ferroelectrics 
(to be published). 

12. K. M. Rabe and U. V. Waghmare, Phil. Trans. R. Soc. Lond. A 354, 2897 (1996). 

13. A. Garcia and D. Vanderbilt, unpublished (cond-mat/9712312). 

14. K. M. Rabe, these Proceedings. 

15. S. Takaoka and K. Murase, Phys. Rev. B 20, 2823 (1979). 

16. D. K. Hohnke, H. Holloway and S. Kaiser, J. Phys. Chem. Solids 33, 2053 (1972). 

17. M. C. Payne et al., "CASTEP 2.1", Cavendish Laboratory, University of Cambridge 
(1991). 

18. U. Waghmare, Ph. D. thesis, Yale University (1996). 

19. If a Te atom belongs to Pb-Te-Ge chain segments in more than one Cartesian direc- 
tion, as for example the z = Te atoms in the oP8 and oC16 configurations, then 
its relaxation is given approximately by the vector sum of the appopriate 0.12 A re- 
laxations for the individual chains. 

20. For technical reasons, the subset of four configurations with cP8 missing could not 
be used in this analysis. The particular values of Ni, N2, and N3 for that set of four 
configurations gives an indeterminate set of linear equations. 

21. The vector representation breaks into irreps according to the point group of the Ge 
site in the given configuration. 

22. Yu. A. Logachev and B. Ya. Moizhes, Sov. Phys. Solid State 19, 1635 (1977). (Fiz. 
Tverd. Tela (Leningrad) 19, 2793 (1977)). 

23. E. Cockayne and K. M. Rabe, unpublished (cond-mat/9712232). 

24. K. M. Rabe and U. V. Waghmare, Phys. Rev. B 52, 13236 (1995). 



